Introduction

This notebook demonstrates how to perform phase and electrochemical assessments starting from a VASP calculation using Python Materials Genomics (pymatgen) and the Materials Project database (via the Materials API). These notebooks are described in detail in

Deng, Z.; Zhu, Z.; Chu, I.-H.; Ong, S. P. Data-Driven First-Principles Methods for the Study and Design of 
Alkali Superionic Conductors. Chem. Mater. 2017, 29 (1), 281–288 DOI: 10.1021/acs.chemmater.6b02648.

If you find these notebooks useful and use the functionality demonstrated, please consider citing the above work.

Let's start by importing some modules and classes that we will be using.

# Uncomment the subsequent lines in this cell to install dependencies for Google Colab.
# !pip install pymatgen==2022.7.19 pymatgen-analysis-diffusion

%matplotlib inline

import json
import re

import matplotlib as mpl
import palettable
from pymatgen.analysis.phase_diagram import (
    CompoundPhaseDiagram,
    PDPlotter,
    PhaseDiagram,
)
from pymatgen.core import Composition, Element
from pymatgen.entries.compatibility import MaterialsProjectCompatibility
from pymatgen.entries.computed_entries import ComputedEntry
from pymatgen.ext.matproj import MPRester
from pymatgen.io.vasp import Vasprun
from pymatgen.util.plotting import pretty_plot

Preparation

We will first read the results from the vasprun.xml output file from our VASP calculations. Only the lowest energy result is used here.

vasprun = Vasprun("aimd_data/vasprun.xml.relax2.gz")
# include structure so proper correction can be applied for oxides and sulfides
entry = vasprun.get_computed_entry(inc_structure=True)

/Users/shyue/repos/pymatgen/pymatgen/io/vasp/inputs.py:1816: UnknownPotcarWarning: POTCAR data with symbol Li_sv does not match any VASP POTCAR known to pymatgen. There is a possibility your POTCAR is corrupted or that the pymatgen database is incomplete.
  warnings.warn(

To construct the phase diagram, we need all entries in the Li-P-S-Cl chemical space. We will use the MPRester class to obtain these entries from the Materials Project via the Materials API.

rester = MPRester()
mp_entries = rester.get_entries_in_chemsys(["Li", "P", "S", "Cl"])

In addition to all the MP entries, here we also load the computed entries of O/S substituted Li-P-O tenary compounds.

with open("aimd_data/lpo_entries.json") as f:
    lpo_data = json.load(f)
lpo_entries = [ComputedEntry.from_dict(d) for d in lpo_data]

Next, we need to combine all the entries and postprocess them using MaterialsProjectCompatibility. This postprocessing step corrects the energies to account for well-known DFT errors, e.g., in the sulfur binding energy.

compatibility = MaterialsProjectCompatibility()
entry = compatibility.process_entry(entry)
entries = compatibility.process_entries([entry] + mp_entries + lpo_entries)

/Users/shyue/repos/pymatgen/pymatgen/entries/compatibility.py:260: UserWarning: sulfide
  warnings.warn(sf_type)

Phase diagram construction

The phase diagram can then be constructed using the PhaseDiagram class, and plotted using the PDPlotter class.

pd = PhaseDiagram(entries)
plotter = PDPlotter(pd)
plotter.show()

We may observe from the above phase diagram that Li6PS5Cl is not a stable phase (red nodes) in the calculated 0K phase diagram.

The pseudo-ternary Li2S-P2S5-LiCl is constructed using the CompoundPhaseDiagram class.

cpd = CompoundPhaseDiagram(
    entries, [Composition("P2S5"), Composition("Li2S"), Composition("LiCl")]
)
cplotter = PDPlotter(cpd, show_unstable=True)
cplotter.show()

Calculating $E_{\rm hull}$ of Li6PS5Cl

We may evaluate the $E_{\rm hull}$ of Li6PS5Cl using the PDAnalyzer.

ehull = pd.get_e_above_hull(entry)
print("The energy above hull of Li6PS5Cl is %.3f eV/atom." % ehull)

The energy above hull of Li6PS5Cl is 0.029 eV/atom.

Electrochemical Stability

The electrochemical stability can be assessed using a similar phase diagram approach, but using the lithium grand potential instead of the internal energy.

First, we need to identify a reference for lithium chemical potential using the bulk Li energy $\mu_{\rm Li}^0$.

li_entries = [e for e in entries if e.composition.reduced_formula == "Li"]
uli0 = min(li_entries, key=lambda e: e.energy_per_atom).energy_per_atom

The PhaseDiagram class provides a quick way to plot the phase diagram at a particular composition (e.g., Li6PS5Cl) as a function of lithium chemical potential called get_element_profile.

el_profile = pd.get_element_profile(Element("Li"), entry.composition)
for i, d in enumerate(el_profile):
    voltage = -(d["chempot"] - uli0)
    print("Voltage: %s V" % voltage)
    print(d["reaction"])
    print("")

Voltage: -0.0 V
4 Li6PS5Cl + 32 Li -> 4 Li3P + 4 LiCl + 20 Li2S

Voltage: 0.8697028629166663 V
4 Li6PS5Cl + 24 Li -> 4 LiP + 4 LiCl + 20 Li2S

Voltage: 0.9324093885416675 V
4 Li6PS5Cl + 21.71 Li -> 0.5714 Li3P7 + 4 LiCl + 20 Li2S

Voltage: 1.1619996604166705 V
4 Li6PS5Cl + 20.57 Li -> 0.5714 LiP7 + 4 LiCl + 20 Li2S

Voltage: 1.2717621691666672 V
4 Li6PS5Cl + 20 Li -> 4 LiCl + 20 Li2S + 4 P

Voltage: 1.7076498754523812 V
4 Li6PS5Cl -> 4 Li3PS4 + 4 LiCl + 4 Li2S

Voltage: 2.1291348952380966 V
4 Li6PS5Cl -> 4 Li3PS4 + LiS4 + 4 LiCl + 7 Li

Voltage: 2.368810398840578 V
4 Li6PS5Cl -> 1.5 LiS4 + 4 LiCl + 2 P2S7 + 18.5 Li

Voltage: 2.886446450666667 V
4 Li6PS5Cl -> 0.5 LiS4 + 2 P2S7 + 4 SCl + 23.5 Li

Voltage: 3.7828703591666635 V
4 Li6PS5Cl -> 2 P2S7 + 2 S + 4 SCl + 24 Li

This element profile can be plotted as a Li evolution versus voltage using matplotlib as follows.

# Some matplotlib settings to improve the look of the plot.
mpl.rcParams["axes.linewidth"] = 3
mpl.rcParams["lines.markeredgewidth"] = 4
mpl.rcParams["lines.linewidth"] = 3
mpl.rcParams["lines.markersize"] = 15
mpl.rcParams["xtick.major.width"] = 3
mpl.rcParams["xtick.major.size"] = 8
mpl.rcParams["xtick.minor.width"] = 3
mpl.rcParams["xtick.minor.size"] = 4
mpl.rcParams["ytick.major.width"] = 3
mpl.rcParams["ytick.major.size"] = 8
mpl.rcParams["ytick.minor.width"] = 3
mpl.rcParams["ytick.minor.size"] = 4


# Plot of Li uptake per formula unit (f.u.) of Li6PS5Cl against voltage vs Li/Li+.

colors = palettable.colorbrewer.qualitative.Set1_9.mpl_colors
plt = pretty_plot(12, 8)

for i, d in enumerate(el_profile):
    v = -(d["chempot"] - uli0)
    if i != 0:
        plt.plot([x2, x2], [y1, d["evolution"] / 4.0], "k", linewidth=3)
    x1 = v
    y1 = d["evolution"] / 4.0
    if i != len(el_profile) - 1:
        x2 = -(el_profile[i + 1]["chempot"] - uli0)
    else:
        x2 = 5.0

    if i in [0, 4, 5, 7]:
        products = [
            re.sub(r"(\d+)", r"$_{\1}$", p.reduced_formula)
            for p in d["reaction"].products
            if p.reduced_formula != "Li"
        ]

        plt.annotate(
            ", ".join(products), xy=(v + 0.05, y1 + 0.3), fontsize=24, color=colors[0]
        )

        plt.plot([x1, x2], [y1, y1], color=colors[0], linewidth=5)
    else:
        plt.plot([x1, x2], [y1, y1], "k", linewidth=3)

plt.xlim((0, 4.0))
plt.ylim((-6, 10))
plt.xlabel("Voltage vs Li/Li$^+$ (V)")
plt.ylabel("Li uptake per f.u.")
plt.tight_layout()